DEA - specific groups of cells
suppressMessages(library(Seurat))
suppressMessages(library(ggplot2))colorClusters <- c("#6fccc1", "#cf1fa7", "#37ee7c", "#da0322", "#ffd94a", "#207ede", "#fea111", "#b09de9",
"#b86fbb", "#150de2", "#ab3f2f", "#d49b6c", "#773aab", "#759474", "#e4b2b8", "#f88a40",
"#74b5cf", "#e6208e", "#0a3535", "#5b37c3", "#cc4025", "#d2dc7c", "#344d76", "#ba8b07",
"#2e10b6", "#4d0cb0", "#e1668a", "#47a58b", "#d7734a", "#ff58c7", "#edceef", "#a21d59",
"#8c8da5", "#ff3e5e", "#688d4b", "#214fca", "#48a8f5", "#752758", "#3b4a4a", "#11674b",
"#9a1bd9", "#16146d", "#7d7277", "#2d8cce", "#2460b8", "#0f9a8f", "#b11e90", "#54ed4f",
"#987217", "#980119", "#cf0606", "#dd182c", "#88a1e4", "#3b89ba", "#12defd", "#feba3e",
"#e31fa5", "#7b3537", "#098a95", "#1f1716", "#df6c94", "#9965ee", "#b438f9", "#0140f5",
"#e6c4a8", "#d94f7c", "#0dae56", "#b86d1d", "#0577ad", "#464551", "#6b0959", "#ccf705",
"#12271d")srt <- readRDS("/mnt/nmorais-nfs/marta/pB_joana/pC_data/srt-all-samples-after-qc.rds")
srt## An object of class Seurat
## 16864 features across 40956 samples within 1 assay
## Active assay: RNA (16864 features, 3000 variable features)
## 3 dimensional reductions calculated: pca, umap, tsnetable(srt$sample)##
## old0 old1 old3 old5 yg0 yg1 yg3 yg5
## 6199 5596 4650 3995 5524 4434 5550 5008sub_srt <- srt[,srt$sample == "old1" | srt$sample == "yg1"]DimPlot(sub_srt, reduction = "tsne", group.by = "RNA_snn_res.0.4", cols = colorClusters) + geom_hline(yintercept= -29, color = "black", size=0.5) +
geom_hline(yintercept= -18, color = "black", size=0.5) +
geom_vline(xintercept= -14, color = "black", size=0.5) +
geom_vline(xintercept= -4, color = "black", size=0.5) +
ggtitle("Eosinophils - group 1")
DimPlot(sub_srt, reduction = "tsne", group.by = "RNA_snn_res.0.4", cols = colorClusters) + geom_hline(yintercept= 9, color = "black", size=0.5) +
geom_hline(yintercept= 15, color = "black", size=0.5) +
geom_vline(xintercept= 3, color = "black", size=0.5) +
geom_vline(xintercept= -2, color = "black", size=0.5) +
ggtitle("Eosinophils - group 2")
srt_eos <- sub_srt[, sub_srt$RNA_snn_res.0.4 == 8]DimPlot(srt_eos, reduction = "tsne", group.by = "RNA_snn_res.0.4", cols = "#b86fbb") +
geom_hline(yintercept= -29, color = "black", size=0.5) +
geom_hline(yintercept= -18, color = "black", size=0.5) +
geom_vline(xintercept= -14, color = "black", size=0.5) +
geom_vline(xintercept= -4, color = "black", size=0.5) +
ggtitle("Eosinophils - group 1")
DimPlot(srt_eos, reduction = "tsne", group.by = "RNA_snn_res.0.4", cols = "#b86fbb") +
geom_hline(yintercept= 9, color = "black", size=0.5) +
geom_hline(yintercept= 15, color = "black", size=0.5) +
geom_vline(xintercept= 3, color = "black", size=0.5) +
geom_vline(xintercept= -2, color = "black", size=0.5) +
ggtitle("Eosinophils - group 2")
xcoordinates <- Embeddings(srt_eos, reduction = "tsne")[,1]
ycoordinates <- Embeddings(srt_eos, reduction = "tsne")[,2]is.group.1 <- (xcoordinates > -14 & xcoordinates < -4) & (ycoordinates > -29 & ycoordinates < -18)
is.group.2 <- (xcoordinates > -2 & xcoordinates < 3) & (ycoordinates > 9 & ycoordinates < 15)Number of cells per group:
table(is.group.1, is.group.2)## is.group.2
## is.group.1 FALSE TRUE
## FALSE 2 172
## TRUE 868 0group <- as.character(srt_eos$RNA_snn_res.0.4)
group[is.group.1] <- "group 1"
group[is.group.2] <- "group 2"
group <- gsub(8, "other", group)srt_eos$group <- group
DimPlot(srt_eos, reduction = "tsne", group.by = "group") 
Idents(srt_eos) <- group
markers <- FindMarkers(srt_eos,
ident.1 = "group 1",
ident.2 = "group 2",
test.use = "MAST")##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!markers <- markers[order(markers$avg_log2FC, decreasing = TRUE),]
markers#saveRDS(markers, "/mnt/nmorais-nfs/marta/pB_joana/pC_data/tables-after-norm/dea-specific-groups-coordinates/eosinophils-subgroup-markers.rds")
#write.table(markers, "/mnt/nmorais-nfs/marta/pB_joana/pC_data/tables-after-norm/dea-specific-groups-coordinates/eosinophils-subgroup-markers.txt", append = FALSE, quote = FALSE)Top violin plots
srt_eos <- srt_eos[,srt_eos$group != "other"]
srt_eos$group <- factor(srt_eos$group, levels = c("group 1", "group 2"))Idents(srt_eos) <- srt_eos$group
VlnPlot(srt_eos, features = rownames(markers)[1:10])
VlnPlot(srt_eos, features = rownames(markers[order(markers$avg_log2FC, decreasing = FALSE),])[1:10])
sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.3.6 SeuratObject_4.1.3 Seurat_4.1.0
##
## loaded via a namespace (and not attached):
## [1] plyr_1.8.7 igraph_1.3.5
## [3] lazyeval_0.2.2 sp_1.5-0
## [5] splines_4.1.2 listenv_0.8.0
## [7] scattermore_0.8 GenomeInfoDb_1.30.1
## [9] digest_0.6.30 htmltools_0.5.3
## [11] fansi_1.0.3 magrittr_2.0.3
## [13] tensor_1.5 cluster_2.1.4
## [15] ROCR_1.0-11 globals_0.16.1
## [17] matrixStats_0.62.0 spatstat.sparse_3.0-0
## [19] prettyunits_1.1.1 colorspace_2.0-3
## [21] ggrepel_0.9.1 xfun_0.34
## [23] dplyr_1.0.10 crayon_1.5.2
## [25] RCurl_1.98-1.9 jsonlite_1.8.3
## [27] progressr_0.11.0 spatstat.data_3.0-0
## [29] survival_3.4-0 zoo_1.8-11
## [31] glue_1.6.2 polyclip_1.10-4
## [33] gtable_0.3.1 zlibbioc_1.40.0
## [35] XVector_0.34.0 leiden_0.4.3
## [37] DelayedArray_0.20.0 future.apply_1.9.1
## [39] SingleCellExperiment_1.16.0 BiocGenerics_0.40.0
## [41] abind_1.4-5 scales_1.2.1
## [43] DBI_1.1.3 spatstat.random_3.0-1
## [45] miniUI_0.1.1.1 Rcpp_1.0.9
## [47] progress_1.2.2 viridisLite_0.4.1
## [49] xtable_1.8-4 reticulate_1.26
## [51] spatstat.core_2.4-4 stats4_4.1.2
## [53] htmlwidgets_1.5.4 httr_1.4.4
## [55] RColorBrewer_1.1-3 ellipsis_0.3.2
## [57] ica_1.0-3 pkgconfig_2.0.3
## [59] farver_2.1.1 sass_0.4.2
## [61] uwot_0.1.14 deldir_1.0-6
## [63] utf8_1.2.2 tidyselect_1.2.0
## [65] labeling_0.4.2 rlang_1.0.6
## [67] reshape2_1.4.4 later_1.3.0
## [69] munsell_0.5.0 tools_4.1.2
## [71] cachem_1.0.6 cli_3.4.1
## [73] generics_0.1.3 ggridges_0.5.4
## [75] evaluate_0.17 stringr_1.4.1
## [77] fastmap_1.1.0 yaml_2.3.6
## [79] goftest_1.2-3 knitr_1.40
## [81] fitdistrplus_1.1-8 purrr_0.3.5
## [83] RANN_2.6.1 pbapply_1.5-0
## [85] future_1.28.0 nlme_3.1-160
## [87] mime_0.12 compiler_4.1.2
## [89] rstudioapi_0.14 plotly_4.10.0
## [91] png_0.1-7 spatstat.utils_3.0-1
## [93] tibble_3.1.8 bslib_0.4.0
## [95] stringi_1.7.8 highr_0.9
## [97] lattice_0.20-45 Matrix_1.5-1
## [99] vctrs_0.5.0 pillar_1.8.1
## [101] lifecycle_1.0.3 spatstat.geom_3.0-3
## [103] lmtest_0.9-40 jquerylib_0.1.4
## [105] RcppAnnoy_0.0.19 data.table_1.14.4
## [107] cowplot_1.1.1 bitops_1.0-7
## [109] irlba_2.3.5.1 GenomicRanges_1.46.1
## [111] httpuv_1.6.6 patchwork_1.1.2
## [113] R6_2.5.1 promises_1.2.0.1
## [115] KernSmooth_2.23-20 gridExtra_2.3
## [117] IRanges_2.28.0 parallelly_1.32.1
## [119] codetools_0.2-18 MASS_7.3-58.1
## [121] assertthat_0.2.1 SummarizedExperiment_1.24.0
## [123] MAST_1.20.0 withr_2.5.0
## [125] sctransform_0.3.5 S4Vectors_0.32.4
## [127] GenomeInfoDbData_1.2.7 hms_1.1.2
## [129] mgcv_1.8-41 parallel_4.1.2
## [131] grid_4.1.2 rpart_4.1.19
## [133] tidyr_1.2.1 rmarkdown_2.17
## [135] MatrixGenerics_1.6.0 Rtsne_0.16
## [137] Biobase_2.54.0 shiny_1.7.2